<!DOCTYPE html>
<html>
<head>
	<meta charset="utf-8"/>
	<link type="text/css" rel="stylesheet" href="shared/css/style.css"/>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "all"} } });
</script>
<script type="text/javascript" src="shared/js/MathJax/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
</head>
<body>

<div style="display:none">
<span class="math">\(\newcommand{\A}{\mat{A}}\)</span>
<span class="math">\(\newcommand{\B}{\mat{B}}\)</span>
<span class="math">\(\newcommand{\C}{\mat{C}}\)</span>
<span class="math">\(\newcommand{\D}{\mat{D}}\)</span>
<span class="math">\(\newcommand{\E}{\mat{E}}\)</span>
<span class="math">\(\newcommand{\F}{\mat{F}}\)</span>
<span class="math">\(\newcommand{\G}{\mat{G}}\)</span>
<span class="math">\(\newcommand{\H}{\mat{H}}\)</span>
<span class="math">\(\newcommand{\I}{\mat{I}}\)</span>
<span class="math">\(\newcommand{\K}{\mat{K}}\)</span>
<span class="math">\(\newcommand{\L}{\mat{L}}\)</span>
<span class="math">\(\newcommand{\M}{\mat{M}}\)</span>
<span class="math">\(\newcommand{\N}{\mat{N}}\)</span>
<span class="math">\(\newcommand{\One}{\mathbf{1}}\)</span>
<span class="math">\(\newcommand{\P}{\mat{P}}\)</span>
<span class="math">\(\newcommand{\Q}{\mat{Q}}\)</span>
<span class="math">\(\newcommand{\Rot}{\mat{R}}\)</span>
<span class="math">\(\newcommand{\R}{\mathbb{R}}\)</span>
<span class="math">\(\newcommand{\S}{\mathcal{S}}\)</span>
<span class="math">\(\newcommand{\T}{\mat{T}}\)</span>
<span class="math">\(\newcommand{\U}{\mat{U}}\)</span>
<span class="math">\(\newcommand{\V}{\mat{V}}\)</span>
<span class="math">\(\newcommand{\W}{\mat{W}}\)</span>
<span class="math">\(\newcommand{\X}{\mat{X}}\)</span>
<span class="math">\(\newcommand{\Y}{\mat{Y}}\)</span>
<span class="math">\(\newcommand{\argmax}{\mathop{\text{argmax}}}\)</span>
<span class="math">\(\newcommand{\argmin}{\mathop{\text{argmin}}}\)</span>
<span class="math">\(\newcommand{\c}{\vec{c}}\)</span>
<span class="math">\(\newcommand{\d}{\vec{d}}\)</span>
<span class="math">\(\newcommand{\e}{\vec{e}}\)</span>
<span class="math">\(\newcommand{\f}{\vec{f}}\)</span>
<span class="math">\(\newcommand{\g}{\vec{g}}\)</span>
<span class="math">\(\newcommand{\mat}[1]{\mathbf{#1}}\)</span>
<span class="math">\(\newcommand{\min}{\mathop{\text{min}}}\)</span>
<span class="math">\(\newcommand{\n}{\vec{n}}\)</span>
<span class="math">\(\newcommand{\p}{\vec{p}}\)</span>
<span class="math">\(\newcommand{\q}{\vec{q}}\)</span>
<span class="math">\(\newcommand{\r}{\vec{r}}\)</span>
<span class="math">\(\newcommand{\transpose}{{\mathsf T}}\)</span>
<span class="math">\(\newcommand{\tr}[1]{\mathop{\text{tr}}{\left(#1\right)}}\)</span>
<span class="math">\(\newcommand{\t}{\vec{t}}\)</span>
<span class="math">\(\newcommand{\u}{\vec{u}}\)</span>
<span class="math">\(\newcommand{\vec}[1]{\mathbf{#1}}\)</span>
<span class="math">\(\newcommand{\x}{\vec{x}}\)</span>
<span class="math">\(\newcommand{\y}{\vec{y}}\)</span>
<span class="math">\(\newcommand{\z}{\vec{z}}\)</span>
<span class="math">\(\renewcommand{\v}{\vec{v}}\)</span>
</div>

<h1 id="geometryprocessing–curvature">Geometry Processing – Curvature</h1>

<blockquote>
<p><strong>To get started:</strong> Fork this repository then issue</p>

<pre><code>git clone --recursive http://github.com/[username]/geometry-processing-curvature.git
</code></pre>
</blockquote>

<h2 id="installationlayoutandcompilation">Installation, Layout, and Compilation</h2>

<p>See
<a href="http://github.com/alecjacobson/geometry-processing-introduction">introduction</a>.</p>

<h2 id="execution">Execution</h2>

<p>Once built, you can execute the assignment from inside the <code>build/</code> by running
on a given mesh:</p>

<pre><code>./curvature [path to mesh.obj]
</code></pre>

<figure>
<img src="images/elephant-curvatures.jpg" alt="Once built you may toggle between showing Gaussian curvature, mean/min/max
curvature and displaying principal curvature directions." />
<figcaption>Once built you may toggle between showing Gaussian curvature, mean/min/max
curvature and displaying principal curvature directions.</figcaption>
</figure>

<h2 id="background">Background</h2>

<p>In this assignment we explore discrete curvature quantities computed on a
surface. These quantities give us local information about a shape. Beyond
inspecting the surface (the extent of this assignment), these quantities become
the building blocks to:</p>

<ul>
<li>define energies to minimize during smoothing/deformation,</li>
<li>identify salient points and curves on the shape, and</li>
<li>provide initial conditions/constraints for <em>remeshing</em>.</li>
</ul>

<p>The fundamental difference between a segment on the real line and a curve is
the introduction of <a href="https://en.wikipedia.org/wiki/Curvature">curvature</a>. This
is quite natural and intuitive. When we draw a 1D object in the plane or in
space we have the freedom to let that object bend. We quantify this &#8220;bending&#8221;
locally as curvature.</p>

<p><a href="https://en.wikipedia.org/wiki/Curvature#Surfaces">Curvature</a> is also the
fundamental difference between a chunk (i.e., subregion) of the <a href="https://en.wikipedia.org/wiki/Plane_(geometry)">Euclidean
Plane</a> and a
<a href="https://en.wikipedia.org/wiki/Surface_(mathematics)">surface</a> that has been
<a href="https://en.wikipedia.org/wiki/Immersion_(mathematics)">immersed</a> in <span class="math">\(\R^3\)</span> (or
elsewhere). Unlike curves, surfaces can bend in each direction at any point.</p>

<p>We start our discussion assuming a smooth surface <span class="math">\(\S\)</span>. We would like to
categorize points on the surface <span class="math">\(\p ∈ \S\)</span> in terms of how the surface bends or
curves locally. </p>

<h3 id="curvatureofplanarcurves">Curvature of planar curves</h3>

<p>Let us briefly recall how
<a href="https://en.wikipedia.org/wiki/Curvature#Precise_definition">curvature</a> is
defined for a <a href="https://en.wikipedia.org/wiki/Plane_curve">planar curve</a>
<span class="math">\(γ:[0,1] → \R²\)</span>.</p>

<p>There are multiple equivalent definitions.</p>

<h4 id="osculatingcircle">Osculating circle</h4>

<p>We can define the <a href="https://en.wikipedia.org/wiki/Tangent">tangent</a> direction at
a point <span class="math">\(\p = γ(s)\)</span> as the limit of the
<a href="https://en.wikipedia.org/wiki/Secant_line">secant</a> formed between <span class="math">\(\p\)</span> and
another point on the curve <span class="math">\(\q = γ(t)\)</span> as <span class="math">\(\q\)</span> approaches <span class="math">\(\p\)</span>:</p>

<figure>
<img src="images/tangent-as-limit-of-secant.gif" alt="" />
</figure>

<p><span class="math">\[
\t(s) = \lim_{\q→\p} \frac{\q-\p}{\|\q-\p\|} = 
\lim_{t→s} \frac{γ(t)-γ(s)}{\|γ(t)-γ(s)\|} = \frac{γ'(s)}{\|γ'(s)\|}.
\]</span>
It always possible, and often convenient, to assume without loss of generality
that <span class="math">\(s\)</span> is an <a href="https://en.wikipedia.org/wiki/Arc_length">arc length
parameterization</a> of the curve <span class="math">\(γ\)</span> so
that <span class="math">\(\|γ'\| = 1\)</span> and therefor the unit tangent vector is simply <span class="math">\(\t(s) =
γ'(s)\)</span>.</p>

<p>In an analogous fashion, we can consider the limit of the
<a href="https://en.wikipedia.org/wiki/Circumscribed_circle">circumcircle</a>
<span class="math">\(C(\q₁,\p,\q₂)\)</span> that passes
through <span class="math">\(\p\)</span> and points <span class="math">\(\q₁\)</span> and <span class="math">\(\q₂\)</span> before and after it on the curve:</p>

<p><span class="math">\[
C(\p) = \lim_{\q₁,\q₂→\p} C(\q₁,\p,\q₂).
\]</span></p>

<figure>
<img src="images/osculating-circle.gif" alt="" />
</figure>

<p>This limit circle is called the <a href="https://en.wikipedia.org/wiki/Osculating_circle">osculating
circle</a> at the point <span class="math">\(\p\)</span> on
the curve <span class="math">\(γ\)</span>. By construction the tangent of the curve and the circle match at
<span class="math">\(\p\)</span>: they&#8217;re both <span class="math">\(γ'\)</span>. The
<a href="https://en.wikipedia.org/wiki/Radius_of_curvature">radius</a> <span class="math">\(R(\p)\)</span> of the
osculating circle <span class="math">\(C(\p)\)</span> at the the point <span class="math">\(\p\)</span> is proportional to how straight
the curve is locally: as the curve becomes more and more straight then the
radius tends toward infinity. This implies that the radius is inversely
proportional to the &#8220;curvy-ness&#8221; of the curve. Hence, the inverse of the radius
is dubbed the curvature:</p>

<p><span class="math">\[
κ(\p) = \frac{1}{R(\p)}.
\]</span>
The radius is a non-negative measure of length with units meters, so the
curvature <span class="math">\(κ\)</span> is an non-negative scalar with units 1/meters. The radius of the
osculating circle can also be written as a limit of the <a href="https://en.wikipedia.org/wiki/Circumscribed_circle#Cartesian_coordinates_from_cross-_and_dot-products">circumcircle
radius</a>:</p>

<p><span class="math">\[
R(\p) = \lim_{\q₁,\q₂→\p} 
  \frac{‖\q₁-\p‖ ‖\p-\q₂‖ ‖\q₂-\q₁‖}
  {2\left| (\q₁-\p) \quad (\p-\q₂)\right|}.
\]</span></p>

<h4 id="signedcurvature">Signed curvature</h4>

<p>Plugging in our arc-length parameterization this reveals that the curvature
(inverse of radius) is equal to the magnitude of change in the tangent or
equivalently the magnitude of second derivative of the curve:</p>

<p><span class="math">\[
κ(s) = \lim_{t→s} \left\| \frac{γ'(t)-γ'(s)}{t-s} \right| = ‖γ''(s)‖.
\]</span></p>

<p>Because we chose the arc-length parameterization, the only change to the
tangent vector <span class="math">\(γ'\)</span> is a change in <em>direction</em> (as opposed to magnitude, since
<span class="math">\(‖γ'‖ := 1\)</span>). This means that the change&#8211;as a vector itself&#8211;is
<em>orthogonal</em> to the tangent. In other words, the change in tangent <span class="math">\(γ''\)</span> points
along the normal direction <span class="math">\(\hat{\n}\)</span>:</p>

<p><span class="math">\[
γ'' ⋅ γ' = 0 \quad → \quad γ'' ⋅ \hat{\n} = ± κ \hat{\n}.
\label{equ:curvature-normal}
\]</span></p>

<p>If we define an orientation to our curve then we can endow the curvature with
a <a href="https://en.wikipedia.org/wiki/Sign_(mathematics)">sign</a> based on whether the
center of the osculating circle lies on the <a href="https://en.wikipedia.org/wiki/Right-hand_rule">left or right
side</a> of the curve. As already
established, the tangent of the osculating circle and the curve agree, so the
vector pointing toward the circle&#8217;s center must be
<a href="https://en.wikipedia.org/wiki/Perpendicular">perpendicular</a> to the tangent:
i.e., in either the positive or negative
<a href="https://en.wikipedia.org/wiki/Normal_(geometry)">normal</a> directions. </p>

<p>If the orientation agrees with increasing the arc-length parameter <span class="math">\(s\)</span>, then the sign can
be
determined by comparing the second derivative vector <span class="math">\(γ''\)</span> to the unit normal
<span class="math">\(\hat{\n} :=
(γ')^⊥\)</span>. The <a href="https://en.wikipedia.org/wiki/Curvature#Signed_curvature"><em><strong>signed
curvature</strong></em></a> at a point <span class="math">\(\p\)</span> is thus given by:</p>

<p><span class="math">\[
\begin{align}
k(\p) &= \text{sign}(γ''(\p)⋅\hat{\n}))\ κ(\p)  \\
      &= γ''(\p) ⋅ \hat{\n}.
\end{align}
\]</span></p>

<h4 id="movingpointanalogy">Moving point analogy</h4>

<p>This definition neatly conforms to our intuition of a curve as the trajectory
of a moving point. Imagine the curved formed by driving along a particular
trajectory <span class="math">\(γ(t)\)</span>, where we really interpret <span class="math">\(t\)</span> as time.</p>

<p>While <span class="math">\(γ'(t)\)</span> corresponds to your velocity vector and <span class="math">\(‖γ'(t)‖\)</span> corresponds to
your speed, the arc-length (re-)parameterization would correspond to having
your friend re-trace your path traveling at a perfectly uniform speed <span class="math">\(‖γ'(s)‖
= 1\)</span>, where your friends &#8220;time&#8221; <span class="math">\(s\)</span> may be different from yours (it may take
longer or shorter depending if you drove fast or slow).</p>

<p>Curvature in the path corresponds to <em>turning</em> and quite literally the amount
by which your friend needs to turn the steering wheel away from the &#8220;straight&#8221;
position: on a straight course, the steering wheel remains at zero-angle
position and the curvature is zero, on a circular course the steering wheel is
fixed at a constant angle in the left or right direction corresponding to
constant positive or negative curvature respectively.</p>

<p>Changing the steering wheel changes the <em>direction</em> of the vehicle&#8217;s velocity.
For your friend driving at constant speed, this is the <em>only</em> change admissible
to the velocity, hence the curvature exactly corresponds to <span class="math">\(γ''(s)\)</span> and to the
steering wheel angle.</p>

<blockquote>
<p>If somebody wants to make a Sega <a href="https://en.wikipedia.org/wiki/Out_Run">Out
Run</a> inspired gif showing a steering
wheel turning next to a little car tracing a curve, I&#8217;ll be very impressed.</p>
</blockquote>

<h4 id="turningnumber">Turning number</h4>

<p>The integrated signed curvature around a <a href="https://en.wikipedia.org/wiki/Curve">closed
curve</a> must be an integer multiple of
<span class="math">\(2π\)</span>:</p>

<p><span class="math">\[
∮ k(s) \ ds = 2π τ,
\]</span>
where τ is an integer called the &#8220;turning number&#8221; of the curve. </p>

<p>This is a bit surprising at first glance. However, in the <em>moving point
analogy</em> a closed curve corresponds to a period trajectory (e.g., driving
around a race-track). When we&#8217;ve made it once around the track, our velocity
direction (e.g., the direction the vehicle is facing) must be pointing in the
original direction. That is, during the course, the car either have turned all the
way around once (<span class="math">\(τ = 1\)</span>) or turned as much clockwise and it did
counter-clockwise (e.g., on a figure 8 course: <span class="math">\(τ=0\)</span>), or made multiple
loops, etc.</p>

<h4 id="discretecurvature">Discrete curvature</h4>

<p>In the discrete world, if a curve is represented as a piecewise-linear chain of
segments, then it&#8217;s natural to associate curvature with vertices: the segments
are flat and therefor contain no curvature.</p>

<p>A natural analog to the definition of curvature as
the derivative of the tangent vector
(i.e., <span class="math">\(k = ‖γ''‖ = ‖\t'‖\)</span>) is to define <em>discrete curvature</em> as the change in
tangent direction between discrete segments meeting at a vertex:</p>

<p><span class="math">\[
k_i =  ∠ (\v_i - (\v_{i-1}-\v_i)) \v_i\v_{i+1} = θ_i,
\]</span>
that is, the signed <a href="https://en.wikipedia.org/wiki/Internal_and_external_angles"><em>exterior
angle</em></a> <span class="math">\(θ_i\)</span> at
the vertex <span class="math">\(\v_i\)</span>.</p>

<figure>
<img src="images/discrete-curvature.svg" alt="" />
</figure>

<p>The turning number theorem for continuous curves finds an <em>immediate</em> analog in
the discrete case. For a closed polygon the discrete signed angles must sum up
to a multiple of <span class="math">\(2π\)</span> in order to close up:</p>

<p><span class="math">\[
∑_{i=1}^n k_i = 2π τ.
\]</span></p>

<p>In this way, we <em>preserve the structure</em> found in the continuous case in our
discrete analog. This structure preservation leads to an understanding of the
exterior angle as an approximation or discrete analog of the <em>locally
integrated</em> curvature.</p>

<p>Alternatively, we could literally fit an circle to the discrete curve based on
local samples and approximate curvature as the inverse radius of the osculating
circle. This curvature measure (in general) will not obey the turning number
theorem, but (conducted properly) it will converge to the pointwise continuous
values under refinement (e.g., as segment length shrinks).</p>

<p>We will explore these two concepts for surfaces, too: discrete analogs that
preserve continuous structures and discretizations that approximate continuous
quantities in the limit.</p>

<h3 id="curvaturesonsurfaces">Curvature(s) on surfaces</h3>

<p>A surface can be curved locally in multiple ways. Consider the difference
between a flat piece of paper, a spherical ping-pong ball and a saddle-shaped
<a href="https://en.wikipedia.org/wiki/Pringles">Pringles chip</a>. The Pringles chip is
the most interesting because it curves &#8220;outward&#8221; in one direction and &#8220;inward&#8221;
in another direction. In this section, we will learn to distinguish and
classify points on a surface based on how it curves in each direction.</p>

<h4 id="normalcurvature">Normal curvature</h4>

<p>The simplest way to extend the curvature that we defined for planar curves to a
surface <span class="math">\(\S\)</span> is to <em>slice</em> the surface through a given point <span class="math">\(\p∈\S\)</span> with a
<a href="https://en.wikipedia.org/wiki/Plane_(geometry)">plane</a> <span class="math">\(\P\)</span> that is parallel
to the <a href="https://en.wikipedia.org/wiki/Normal_(geometry)">surface normal</a>
<span class="math">\(\n(\p)\)</span>.</p>

<p>The (local) intersection of the surface <span class="math">\(\S\)</span> and the plane <span class="math">\(\P\)</span> will trace a
curve <span class="math">\(γ\)</span>, upon which we can immediately use the planar curvature definition
above. </p>

<figure>
<img src="images/normal-curvature.svg" alt="Source" />
<figcaption><a href="http://brickisland.net/cs177fa12/?p=214">Source</a></figcaption>
</figure>

<p>There are infinitely many planes that pass through a given point <span class="math">\(\p\)</span> and lie
parallel to a given normal vector <span class="math">\(\n(\p)\)</span>: the plane can rotate around the
normal <span class="math">\(\n(\p)\)</span> by any angle <span class="math">\(φ\)</span>. For each choice of <span class="math">\(φ\)</span>, the plane will define
an intersecting curve <span class="math">\(γ_φ\)</span> and thus for every angle <span class="math">\(φ\)</span> there will be a
<a href="https://en.wikipedia.org/wiki/Curvature#Normal_sections"><em>normal curvature</em></a>:</p>

<p><span class="math">\[
k_\n(φ,\p) = γ''_\n(\p).
\]</span></p>

<h4 id="meancurvature">Mean curvature</h4>

<p>Normal curvature requires choosing an angle, so it doesn&#8217;t satiate our
desire to reduce the &#8220;curvy-ness&#8221; to a single number for any point on the
surface. A simple way to reduce this space of normal curvatures is to, well,
average all possible normal curvatures. This defines the <a href="https://en.wikipedia.org/wiki/Mean_curvature">mean
curvature</a>:</p>

<p><span class="math">\[
H(\p) = \frac{1}{2π}∫_0^{2π} k_\n(φ,\p) \ dφ.
\]</span></p>

<h4 id="maximumandminimumcurvature">Maximum and minimum curvature</h4>

<p>Another obvious way to reduce the space of normal curvatures to a single number
is to consider the maximum or minimum normal curvature over all choices of <span class="math">\(φ\)</span>:</p>

<p><span class="math">\[
\begin{align}
k₁(\p) &= \max_φ \ k_\n(φ,\p) \\
k₂(\p) &= \min_φ \ k_\n(φ,\p).
\end{align}
\]</span></p>

<p>Collectively, these are referred to as the <a href="https://en.wikipedia.org/wiki/Principal_curvature">principal
curvatures</a> and
correspondingly the angles that maximize and minimize curvature are referred to
as the principle curvature directions:</p>

<p><span class="math">\[
\begin{align}
φ₁(\p) &= \argmax_φ \ k_\n(φ,\p) \\
φ₂(\p) &= \argmin_φ \ k_\n(φ,\p).
\end{align}
\]</span></p>

<p><a href="https://en.wikipedia.org/wiki/Euler%27s_theorem_(differential_geometry)">Euler&#8217;s
theorem</a>
states that the normal curvature is a quite simple function of <span class="math">\(φ\)</span> and the
principle curvatures:</p>

<p><span class="math">\[
k_\n(φ,\p) = k₁ \cos^2 φ + k₂ \sin^2 φ.
\]</span></p>

<p>There are two immediate and important consequences:</p>

<ol>
<li>the principal curvature directions (<span class="math">\(φ₁\)</span> and <span class="math">\(φ₂\)</span>) are orthogonal, and</li>
<li>the mean curvature reduces to the average of principal curvatures:</li>
</ol>

<p><span class="math">\[
H = ½(k₁ + k₂).
\]</span></p>

<blockquote>
<p>For more theory and a proof of Euler&#8217;s theorem, I recommend &#8220;Elementary
Differential Geometry&#8221; by Barret O&#8217;Neill, Chapter 5.2.</p>
</blockquote>

<h4 id="gaussiancurvature">Gaussian curvature</h4>

<p>Maximum, minimum and mean curvature reduce curvature to a single number, but
still cannot (alone) distinguish between points lying on a round ping-pong
ball, a flat sheet of paper, the cylindrical Pringles can and a saddle-shaped
Pringles chip.</p>

<p>The neck of this cartoon elephant&#8211;like a Pringles chip&#8211;bends inward in one
direction (positive <span class="math">\(k₁ &gt; 0\)</span>) and outward in the other
direction (negative <span class="math">\(k₂ &lt; 0\)</span>). </p>

<figure>
<img src="images/cartoon-elephant-principal-and-gaussian-curvature.jpg" alt="Maximum $k₁$, minimum $k₂$,and Gaussian curvature $K = k₁k₂$" />
<figcaption>Maximum <span class="math">\(k₁\)</span>, minimum <span class="math">\(k₂\)</span>,and Gaussian curvature <span class="math">\(K = k₁k₂\)</span></figcaption>
</figure>

<p>The <em>product</em> of the principal curvatures maintains the disagreement in sign
that categories this saddle-like behavior. This product is called <a href="https://en.wikipedia.org/wiki/Gaussian_curvature">Gaussian
curvature</a>:</p>

<p><span class="math">\[
K = k_1 k_2.
\]</span></p>

<h4 id="relationshiptosurfacearea">Relationship to surface area</h4>

<p>Both mean and Gaussian curvature have meaningful relationships to surface
area.</p>

<h5 id="meancurvatureasareagradient">Mean Curvature as area gradient</h5>

<p>Let us consider a seemingly unrelated yet familiar problem. Suppose we would
like to minimize the surface area of a given shape <span class="math">\(\S\)</span>. </p>

<p>If we consider a small patch <span class="math">\(A ⊂ \S\)</span> centered around some point <span class="math">\(\p ∈\S\)</span>, then
we can imagine forces at the boundary of the patch <span class="math">\(∂A\)</span> that try to shrink the
area of the patch. These &#8220;surface tension&#8221; forces point inward against the
normal direction of the patch boundary <span class="math">\(η\)</span>. Consider the integral effect of
these forces and apply <a href="https://en.wikipedia.org/wiki/Divergence_theorem">divergence
theorem</a> to move the integral
to the interior of the patch.</p>

<p><span class="math">\[
∮_{∂A} -η \;ds = ∮_{∂A} \n × d\x = ∫_A ∆\x \;dA,
\]</span>
where <span class="math">\(\n\)</span> is the surface normal, <span class="math">\(d\x\)</span> is the infinitesimal change in position
around the boundary of the patch, and <span class="math">\(∆\x ∈ \R³\)</span> computes the
<a href="https://en.wikipedia.org/wiki/Laplace_operator">Laplacian</a> of the position
coordinate functions.</p>

<p>In the limit as <span class="math">\(A\)</span> shrinks around any point on the surface,
this relationship gives us the intuition that <em>moving</em> the surface in the
direction <span class="math">\(∆\x\)</span> will decrease the surface area. </p>

<p>The Laplacian <span class="math">\(∆f\)</span> of a function <span class="math">\(f\)</span> on the surface does not depend on the
choice of parameterization. It is defined as the divergence of the gradient of
the function or equivalently the trace of the Hessian:</p>

<p><span class="math">\[
∆f = ∇⋅ ∇f = \tr{ 
\left[ 
\begin{array}{cc}
\frac{∂²f}{∂u²} & \frac{∂²f}{∂u∂v} \\
\frac{∂²f}{∂v∂u} & \frac{∂²f}{∂v²} 
\end{array}
\right]
}
=
\frac{∂²f}{∂u²} + \frac{∂²f}{∂v²}.
\]</span></p>

<p>If we generously choose <span class="math">\(u\)</span> and <span class="math">\(v\)</span> to vary in the principal directions <span class="math">\(φ₁\)</span>
and <span class="math">\(φ₂\)</span> above. In this case, the Laplacian <span class="math">\(∆\x\)</span> of the position function
reduces to the sum of principal curvatures times the normal (recall Equation
<span class="math">\((\ref{equ:curvature-normal})\)</span>):</p>

<p><span class="math">\[
\begin{align}
∆\x &= \frac{∂²\x}{∂u²} + \frac{∂²\x}{∂v²} \\
    &= k₁ \n + k₂ \n \\
    &= 2H\n,
\end{align}
\]</span>
where <span class="math">\(H\n ∈ \R³\)</span> is called the <em><strong>mean curvature normal</strong></em> vector. We have
shown that the mean curvature normal is equal half the Laplacian of the
embedding function.</p>
<!--
This turns out to be surprisingly complicated to derive without "going of the
deep end" into differential forms etc.

The trivial definition of surface area as an integral over $\S$
becomes an integral over $Ω$ of the determinant of the
[Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant) of
$\x$ after applying [substitution
rule](https://en.wikipedia.org/wiki/Integration_by_substitution):

\\[
∫_\S 1\ dA = ∫_Ω \left| ∇\x \right| dΩ
\\]

-->

<h5 id="gaussiancurvatureasareadistortion">Gaussian Curvature as area distortion</h5>

<p>As the product of principal curvatures, Gaussian curvature <span class="math">\(K = k₁k₂\)</span> measures zero
anytime one (or both) of the principal curvatures are zero. Intuitively, this
happens only for surfaces that curve or bend in one direction. Imagine rolling
up a sheet of paper. Surfaces with zero Gaussian curvature <span class="math">\(K = 0\)</span> are called
<em>developable surfaces</em> because the can be flattened (developed) on to the flat
plane (just as you might unroll the piece of paper) <em>without</em> stretching or
shearing. As a corollary, surfaces with non-zero Gaussian curvature <em>cannot</em> be
flattened to the plane without stretching some part.</p>

<p>Locally, Gaussian curvature measures how far from developable the surface is:
how much would the local area need to stretch to become flat.</p>

<p>First, we introduce the <a href="https://en.wikipedia.org/wiki/Gauss_map">Gauss map</a>, a
continuous map <span class="math">\(N:\S→S²\)</span> from every point <span class="math">\(\p\)</span> on the surface <span class="math">\(\S\)</span> to the unit
sphere <span class="math">\(S²\)</span> so that <span class="math">\(N(\p) := \n(\p)\)</span>, the unit normal at <span class="math">\(\p\)</span>.</p>

<p>Consider a small patch on a curved surface. Gaussian curvature <span class="math">\(K\)</span> can
equivalently be defined as the limit of the ratio between the area
area <em>swept</em> out by the unit normal on the Gauss map <span class="math">\(A_G\)</span> and
the area of the surface patch <span class="math">\(A\)</span>:</p>

<p><span class="math">\[
K = \lim_{A→0} \frac{A_G}{A}.
\label{equ:gaussian-curvature-area}
\]</span></p>

<p>Let&#8217;s consider different types of regions:</p>

<ul>
<li>flat: <span class="math">\(A_G=0\)</span> because the Gauss map is a point,</li>
<li>cylindrical: <span class="math">\(A_G=0\)</span> because the Gauss map is a curve,</li>
<li>spherical: <span class="math">\(A_G&gt;0\)</span> because the Gauss map will maintain positive swept-area,
 and</li>
<li>saddle-shaped: <span class="math">\(A_G&lt;0\)</span> because the area on the Gauss map will maintain
 <em>oppositely</em> oriented area (i.e., from the spherical case).</li>
</ul>

<p><img src="images/gauss-map-plane.gif" alt="A patch on a plane and its corresponding patch on the Gauss map." />
<img src="images/gauss-map-cylinder.gif" alt="A patch on a cylinder and its corresponding patch on the Gauss map." />
<img src="images/gauss-map-sphere.gif" alt="A patch on a sphere and its corresponding patch on the Gauss map." />
<img src="images/gauss-map-saddle.gif" alt="A patch on a saddle and its corresponding patch on the Gauss map." /></p>

<p>Similar to the turning number theorem for curves, there exists an analogous
<a href="https://en.wikipedia.org/wiki/Gauss–Bonnet_theorem">theorem for surfaces</a>
stating that the total Gaussian curvature must be an integer multiple of <span class="math">\(2π\)</span>:</p>

<p><span class="math">\[
∫_S K dA = 2π χ(\S),
\label{equ:gauss-bonnet}
\]</span>
where <span class="math">\(χ(\S)\)</span> is the <a href="https://en.wikipedia.org/wiki/Euler_characteristic">Euler
characteristic</a> of the
surfaces <span class="math">\(\S\)</span> (a topological <em>invariant</em> of the surface revealing how many
<a href="https://en.wikipedia.org/wiki/Genus_(mathematics)">holes</a> the surface has).</p>

<p>In stark contrast to mean curvature, this theorem tells us that we cannot add
Gaussian curvature to a surface without:</p>

<ol>
<li>removing an equal amount some place else, or</li>
<li>changing the topology of the surface.</li>
</ol>

<p>Since changing the topology of the surface would require a discontinuous
deformation, adding and removing Gaussian curvature must also balance out for
smooth deformations. This simultaneously explains why a cloth must have
wrinkles when draping over a table, and why a deflated basketball will not lie
flat on the ground.</p>

<h4 id="shapeoperator">Shape operator</h4>

<p>There is yet another way to arrive at principal, mean and Gaussian curvatures.
Consider a point <span class="math">\(\p\)</span> on a surface <span class="math">\(\S\)</span> with unit normal vector <span class="math">\(\n\)</span>. If we
pick a unit tangent vector <span class="math">\(\v\)</span> (i.e., so that <span class="math">\(\v ⋅ \n = 0\)</span>), then we can ask
how does the normal <span class="math">\(\n\)</span> change as we move in the direction of <span class="math">\(\v\)</span> along the
surface:</p>

<p><span class="math">\[
S_\p(\v) := ∇ \n ⋅ \v
\]</span>
we call <span class="math">\(S_\p\)</span> the <a href="https://en.wikipedia.org/wiki/Differential_geometry_of_surfaces#Shape_operator"><em><strong>shape
operator</strong></em></a>
at the point <span class="math">\(\p\)</span>. Just as how in Equation <span class="math">\((\ref{equ:curvature-normal})\)</span> the
curvature normal must point in the normal direction, the shape operator takes
as input a tangent vector and outputs another tangent vector (i.e., the change
in the unit normal must be tangential to the surface; no change can occur in
the normal direction itself).</p>

<p>Locally, the tangent vector space is two-dimensional spanned by basis vectors
<span class="math">\(\e₁,\e₂ ∈ \R²\)</span> so we can think of the
shape operator as a mapping from <span class="math">\(\R²\)</span> to <span class="math">\(\R²\)</span>. As a differential operator,
the shape operator is a <em>linear operator</em>. This means we can represent its
action on a tangent vector <span class="math">\(\v = x₁ \e₁ + x₂\e₂\)</span> as a matrix:</p>

<p><span class="math">\[
S_\p(\v) = 
\left[
\begin{array}{cc}
S_\p(\e₁)⋅\e₁ & S_\p(\e₁)⋅\e₂ \\
S_\p(\e₂)⋅\e₁ & S_\p(\e₂)⋅\e₂
\end{array}
\right] \v
\]</span></p>

<p>Given <span class="math">\(\r₁\)</span> and <span class="math">\(\r₂\)</span> are the principal curvature directions (as unit 2D tangent
vectors) we can rotate our coordinate frame to align <span class="math">\(\e₁\)</span> and <span class="math">\(\e₂\)</span> with the
principal curvature directions. The shape operator takes on a very special
form:</p>

<p><span class="math">\[
S_\p = 
\left[\r₁ \quad \r₂\right]
\left[
\begin{array}{cc}
k₁ & 0 \\
0 & k²
\end{array}
\right]
\left[\r₁ \quad \r₂\right]^\transpose
\]</span></p>

<blockquote>
<p>Consider why the off-diagonal terms are zero. Think about the <em>extremality</em>
of the principal curvatures.</p>
</blockquote>

<p>We have actually conducted an <a href="https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">eigen
decomposition</a> on
the shape operator. Reading this progression backwards, the eigen decomposition
of the shape operator expressed in any basis will reveal:</p>

<ol>
<li>the principal curvatures as the eigen values, and</li>
<li>the principal curvature directions as the eigen vectors.</li>
</ol>
<!--
Curvature for curves is the change in tangent vector under an arc-length. For a
given tangent direction on a surface, we extended this definition to define the
normal curvature as the curvature of the curve made by interesting surface with
a plane aligned with the chosen direction. Given two orthonormal tangent
directions $\u$ and $\v$ (i.e., a local parameterization), let's collect the
normal curvature normal vectors:

\\[
\begin{align}
k_\n(\u,\p)\n &= γ''_\u(\p) \\
k_\n(\v,\p)\n &= γ''_\v(\p).
\end{align}
\\]

If instead we equivalent consider the _change in normal vector_ for each sliced
curve, our curvature vectors will live in the orthogonal space: the tangent
space. 

\\[
\begin{align}
k_\n(\u,\p)\v &= γ''_\u(\p) \\
k_\n(\v,\p)\u &= γ''_\v(\p).
\end{align}
\\]

> Before we chose the normal direction by an angle $φ$, but for any tangent
> direction $\u$ we can determine its correspond $φ$ so that $k_\n(φ) =
> k_\n(\u)$.
-->

<h3 id="discretecurvaturesonsurfaces">Discrete curvatures on surfaces</h3>

<h4 id="discretemeancurvaturenormalviadiscretelaplace">Discrete mean curvature normal via discrete Laplace</h4>

<p>By now we are very familiar with the discrete Laplacian for triangle meshes:</p>

<p><span class="math">\[
∆f ≈ \M^{-1} \L \f,
\]</span>
where <span class="math">\(\M,\L ∈\R^{n×n}\)</span> are the mass and cotangent matrices respectively.</p>

<p>When applied to the vertex positions, this operator gives a point-wise (or
rather integral average) approximation of the mean curvature normal:</p>

<p><span class="math">\[
H\n ≈ \H = \M^{-1} \L \V ∈ \R^{n×3}.
\]</span></p>

<p>Stripping the magnitude off the rows of the resulting matrix would give the
<em>unsigned</em> mean curvature. To make sure that the sign is preserved we can check
whether each row in <span class="math">\(\H\)</span> agrees or disagrees with consistently oriented
per-vertex normals in <span class="math">\(\N ∈ \R^{n×3}\)</span>. </p>

<p>This connection between the Laplace operator and the mean curvature normal
provides additional understanding for its use as a geometric smoothing operator
(see &#8220;Computing Discrete Minimal Surfaces and Their Conjugates&#8221; [Pinkall and
Polthier 1993]).</p>

<h4 id="discretegaussiancurvatureviaangledefect">Discrete Gaussian curvature via angle defect</h4>

<p>On a discrete surface represented as a triangle mesh, curvature certainly can&#8217;t
live on the flat faces. Moreover, Gaussian curvature can&#8217;t live along edges
because we can always <em>develop</em> the triangles on either side of an edge to the
plane without stretching them. In fact we can develop any arbitrarily long
chain of faces connected by edges so long as it doesn&#8217;t form a loop or contain
all faces incident on a vertex. This hints that discrete Gaussian curvature
(like curvature for curves) must live at vertices.</p>

<p>Using the definition of Gaussian curvature in terms of the area on the Gauss
map in Equation <span class="math">\((\ref{equ:gaussian-curvature-area})\)</span>: flat faces correspond
points on the Gauss map (contributing nothing), edges correspond to area-less
curves (traced by their <a href="https://en.wikipedia.org/wiki/Dihedral_angle">dihedral
angles</a>), but vertices correspond
to spherical polygons connecting face normal-points. The area <span class="math">\(Ω\)</span> subtended on
the Gauss map is call the <a href="https://en.wikipedia.org/wiki/Solid_angle">solid
angle</a>. Conveniently, this area is
simply the <a href="https://en.wikipedia.org/wiki/Angular_defect#Descartes.27_theorem">angle
defect</a> of
internal angles <span class="math">\(θ_f\)</span> incident on the <span class="math">\(i\)</span>-th vertex contributed by each <span class="math">\(f\)</span>-th
incident face:</p>

<p><span class="math">\[
Ω_i = 2π - ∑\limits_{f ∈ \text{faces(i)}} θ_{if}.
\]</span></p>

<figure>
<img src="images/angle-defect.png" alt="Gaussian Curvature and Shell Structures " />
<figcaption>&#8220;Gaussian Curvature and Shell Structures&#8221; [Calladine
1986]</figcaption>
</figure>

<p>Thus, our discrete analog of locally <em>integrated</em> Gaussian curvature is given
as the angle defect at the <span class="math">\(i\)</span>-th vertex. The local integral average (or
<em>pointwise</em>) discrete Gaussian curvature is the angle defect divided by the
local area associated with the <span class="math">\(i\)</span>-th vertex:</p>

<p><span class="math">\[
K_i = \frac{2π - ∑\limits_{f ∈ \text{faces(i)}} θ_{if}}{A_i}.
\]</span></p>

<p>By way of closing up the Gauss map, closed polyhedral surfaces (i.e., meshes)
will obey the
<a href="https://en.wikipedia.org/wiki/Gauss–Bonnet_theorem">Gauss-Bonnet</a> in Equation
<span class="math">\((\ref{equ:gauss-bonnet})\)</span>, too:</p>

<p><span class="math">\[
∑\limits_{i=1}^n K_i = 2π χ(\S).
\]</span></p>

<p>We can connect this to <a href="https://en.wikipedia.org/wiki/Euler_characteristic">Euler&#8217;s
formula</a> for polyhedra in our very first
assignment:</p>

<p><span class="math">\[
\frac{1}{2π} ∑\limits_{i=1}^n K_i =  |V| - |E| + |F|,
\]</span>
where <span class="math">\(|V|, |E|, |F|\)</span> are the number of vertices, edges and faces respectively.</p>

<h4 id="approximationandeigendecompositionoftheshapeoperator">Approximation and eigen decomposition of the shape operator</h4>

<p>Alternatively, we can approximate all curvatures of a surface by locally
fitting an analytic surface and <em>reading</em> off its curvature values. Since
planes have no curvature, the simplest type of analytic surface that will give
a non-trivial curvature value is a quadratic surface.</p>

<p>Thus, the algorithm proceeds as follows. For each vertex <span class="math">\(\v\)</span> of the given mesh,</p>

<ol>
<li>gather a sampling of points in the vicinity. For simplicity, let&#8217;s just
 grab all other vertices that share an edge with <span class="math">\(\v\)</span> or share an edge with a
 vertex that shares an edge with <span class="math">\(\v\)</span> (i.e., the &#8220;two-ring&#8221; of <span class="math">\(\v\)</span>). For most
 sane meshes, this will provide enough points. Gather the positions of these
 <span class="math">\(k\)</span> points <em>relative</em> to <span class="math">\(\v\)</span> (i.e., <span class="math">\(\v_i - \v\)</span>) into a matrix <span class="math">\(\P ∈
 \R^{k×3}\)</span>.</li>
<li>Next, we are going to define a quadratic surface as a height field above
 some two-dimensional plane passing through <span class="math">\(\v\)</span>. Ideally, the plane is
 orthogonal to the normal at <span class="math">\(\v\)</span>. To find such a plane, compute the
 <a href="https://en.wikipedia.org/wiki/Principal_component_analysis">principal-component
 analysis</a> of <span class="math">\(\P\)</span>
 (i.e., conduct eigen decomposition on <span class="math">\(\P^\transpose \P\)</span>). Let <span class="math">\(\S ∈ \R^{k ×
 2}\)</span> be the coefficients for two most principal directions (call them the <span class="math">\(u\)</span>-
 and <span class="math">\(v\)</span>- directions) corresponding to each point in <span class="math">\(\P\)</span>, and let <span class="math">\(\B ∈
 \R^{k}\)</span> be the &#8220;height&#8221; of each point in the least principal direction (call
 it the <span class="math">\(w\)</span>-direction).</li>
<li>An quadratic function as a height-field surface passing through the origin
 is given by:
<span class="math">\[
w = a₁ u + a₂ v + a₃ u² + a₄uv + a₅v².
\]</span>
 We have <span class="math">\(k\)</span> sets of <span class="math">\(u,v\)</span> values and <span class="math">\(w\)</span> values. Treat this as a
 least-squares fitting problem and solve for the 5 unknown coefficients.
 (<code>igl::pinv</code> is good for solving this robustly).</li>
<li>Each element of the shape operator for the graph of a quadratic function
 over the plane has a closed form expression. You need to derive these by hand.
 Just kidding. The shape operator can be constructed as the product of two
 matrices:
<span class="math">\[
S = -
\left[
\begin{array}{cc}
e & f \\
f & g
\end{array}
\right]
\left[
\begin{array}{cc}
E & F \\
F & G
\end{array}
\right]^{-1}
\]</span>
 known as the second and first fundamental forms respectively. The entries of
 these matrices categorize the stretch and bending in each direction:
<span class="math">\[
E = 1+a₁² \\
F = a₁a₂  \\
G = 1+a₂² \\
e = \frac{2a₃}{\sqrt{a₁² + 1 + a₂²}} \\
f = \frac{a₄}{\sqrt{a₁² + 1 + a₂²}} \\
g = \frac{2a₅}{\sqrt{a₁² + 1 + a₂²}}
\]</span>
 See Table 1 of &#8220;Estimating Differential Quantities Using Polynomial Fitting of
 Osculating Jets&#8221; [Cazals &amp; Pouget 2003] to double check for typos :-).</li>
<li>Eigen decomposition of <span class="math">\(S\)</span> reveals the principal curvatures <span class="math">\(k₁\)</span> and <span class="math">\(k₂\)</span>
 <em>and</em> the principal tangent directions (in the <span class="math">\(uv\)</span> PCA basis).</li>
<li>Lift the principal tangent directions back to world <span class="math">\(\R³\)</span> coordinates.</li>
</ol>

<h2 id="tasks">Tasks</h2>

<p><a href="https://archive.org/details/ElementaryDifferentialGeometry">Download</a> Barret
O&#8217;Neill&#8217;s book. This is my go-to differential geometry book. The section on
curvature and the shape operator should help resolve questions and fill in
missing proofs above.</p>

<h3 id="blacklist">Blacklist</h3>

<ul>
<li><code>igl::gaussian_curvature</code></li>
<li><code>igl::internal_angles</code> (or any of the other overloads)</li>
<li><code>igl::principal_curvatures</code></li>
</ul>

<h3 id="whitelist">Whitelist</h3>

<ul>
<li><code>igl::adjacency_matrix.h</code></li>
<li><code>igl::cotmatrix</code></li>
<li><code>igl::invert_diag</code></li>
<li><code>igl::massmatrix</code></li>
<li><code>igl::per_vertex_normals</code></li>
<li><code>igl::pinv</code></li>
<li><code>igl::slice</code></li>
<li><code>igl::sort</code></li>
<li><code>igl::squared_edge_lengths</code></li>
</ul>

<h3 id="srcmean_curvature.cpp"><code>src/mean_curvature.cpp</code></h3>

<p>Compute the discrete mean curvature at each vertex of a mesh (<code>V</code>,<code>F</code>) by
taking the signed magnitude of the mean curvature normal as a <em>pointwise</em> (or
<em>integral average</em>) quantity.</p>

<h3 id="srcinternal_angles.cpp"><code>src/internal_angles.cpp</code></h3>

<p>Given (squared) edge-lengths of a triangle mesh <code>l_sqr</code> compute the internal
angles at each corner (a.k.a. wedge) of the mesh.</p>

<h3 id="srcangle_defect.cpp"><code>src/angle_defect.cpp</code></h3>

<p>Compute the discrete angle defect at each vertex of a triangle mesh
(<code>V</code>,<code>F</code>), that is, the <em>locally integrated</em> discrete Gaussian
curvature.</p>

<h3 id="srcprincipal_curvatures.cpp"><code>src/principal_curvatures.cpp</code></h3>

<p>Approximate principal curvature values and directions locally by considering
the two-ring neighborhood of each vertex in the mesh (<code>V</code>,<code>F</code>).</p>

</body>
</html>
